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Abstract. We show how the Gillespie algorithm, originally developed to describe 
coupled chemical reactions, can be used to perform numerical simulations of a granular 
intruder particle colliding with thermalized bath particles. The algorithm generates 
a sequence of collision "events" separated by variable time intervals. As input, it 
requires the position-dependent flux of bath particles at each point on the surface of 
the intruder particle. We validate the method by applying it to a one-dimensional 
system for which the exact solution of the homogeneous Boltzmann equation is known 
and investigate the case where the bath particle velocity distribution has algebraic tails. 
We also present an application to a granular needle in bath of point particles where 
we demonstrate the presence of correlations between the translational and rotational 
degrees of freedom of the intruder particle. The relationship between the Gillespie 
algorithm and the commonly used Direct Simulation Monte Carlo (DSMC) method is 
also discussed. 
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1. Introduction 

Kinetic theories of granular systems are usually constructed starting from the Boltzmann 
equation or one of its variants[l[ 0|. Rarely, it is possible to obtain an exact, analytic 
solution of the Boltzmann equation^. More typically, however, approximations are 
required. It is then highly desirable to assess the quality of the theoretical prediction 
by comparing it with accurate numerical solutions of the Boltzmann equation. It is 
the purpose of this paper to show that, besides the celebrated Direct Simulation Monte 



Carlo (DSMC) introduced by Bird 4], there exists an alternative method, originally 
proposed by Gillespie ^; to study coupled chemical reactions. 

One class of system that is amenable to a Boltzmann approach and that has received 
considerable attention in recentyears consists of a single intruder (or tracer) particle in 
a bath of thermalized particles [7; BSE HQ, 

showing in particular the absence of 
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equipartition. This phenomena has been observed experimentally in two- dimensional jl3| 



and three-dimensional 14J granular gases. The intruder-bath particle collisions are 
dissipative, while the bath particles have an ideal gas structure and a specified velocity 
distribution characterized by a fixed temperature. Since Gaussian velocity statistics are 
quite rare, it is important that numerical and theoretical approaches be able to treat a 
general distribution. 

We begin by outlining the Gillespie algorithm and how it can be used to obtain 
a numerical solution of the Boltzmann equation. We then illustrate the application 
of the algorithm with two examples. The first is a one- dimensional system consisting 
of an intruder particle in a bath of thermalized point particles. If the bath particles 
have a Gaussian velocity distribution one recovers the exact Gaussian intruder velocity 
distribution function p with a granular temperature that is smaller than the bath 
temperature. In addition when we impose a power-law velocity distribution on the bath 
particles, we find the same form for the intruder particle distribution. 

The second system is two-dimensional and consists of a needle intruder in a bath of 
point particles. It is known that equipartition does not hold between different degrees of 
freedom of the needle [lfl| . Here we use the Gillespie algorithm to obtain a new physical 
result, namely the presence of correlations between the translational and rotational 
degree of freedom of the needle. 

2. Algorithm 

The model consists of a single intruder particle that undergoes a series of collisions with 
the surrounding bath particles. Let P(t) denote the probability that no event (collision) 
has occurred in the interval (0, t). Then 

Pit + At) = P(t){l - <f>(t)At + 0(At 2 )) (1) 

where <p(t) is the event rate in general and the collision rate in this application. 
Expanding the lhs to first order in At and taking the limit At — > leads to 

T = < 2) 

Integrating and using the boundary condition that P(0) = 1 gives 



P(t) = exp[- / (j){t')dt'\ (3) 
Jo 

If <p(t) is constant between collisions this expression takes the simple form: 

P(t) = exp(-0t) (4) 

At the end of the waiting time an event (collision) occurs that alters the value of <p or 
the way that evolves with time. We now consider two specific applications. In the 
first, the flux of colliding particles is constant between collisions and the simpler form, 
Eq |U may be used. Our second example, an anisotropic object that rotates between 
collisions, requires the use of the more general form, Eq El 
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3. Applications 

3.1. One- dimensional system 

Consider first a one-dimensional system consisting of an intruder particle of mass M 
moving in a bath of thermalized point particles each of mass m. The dynamics of the 
intruder is described by the Boltzmann equation. 

The velocity distribution of the bath particles is denoted by f(v,a), where a is 
related to the bath temperature, Tb- 

The flux of particles that collide with the right hand side of the intruder particle 
moving with a velocity V\ is: 

/VI 
(v!-v)f(v,a)dv, (5) 
-oo 

where p is the number density of the bath particles. Similarly the flux on the left hand 
side is: 

POO 

0-Oi)=P/ (v -v 1 )f(v,a)dv (6) 
and the total collision rate is 

( f>(v 1 )=ct> + (vi) + <t>-(vi) (7) 

Let F(vi,t) denote the time- dependent distribution function of tracer particle 
velocity. Then this evolves according to 

= -<j>{vi)F{v u t) + J°° ij(v -> Vl )F(v, t)dv (8) 

which is equivalent to the homogeneous Boltzmann equation^. The first term on the 
rhs is a loss term corresponding to the probability that a tracer particle with velocity 
v i undergoes a collision (necessarily to a different velocity) per unit time. The second, 
or gain, term contains the function ip(v — > Vi) that is the rate that tracer particles 
moving with velocity v are transformed (by collisions with the bath particles) to those 
moving with velocity V\ . For collisions with the rhs of the intruder particle the explicit 
expression is 

, . , [I + mV, , „, 1+M, , \ 

ip+(v -»• v x ) = p — — {v-v 1 )f(v + — (ui-<u),a), (9) 

\ 1 + a J 1 + a 

with vi < v. A similar expression applies for collisions on the left hand side of the 
intruder particle. We note that a representation of the Boltzmann equation similar to 
EqlHlwas employed by Puglisi et al. 

The Gillespie algorithm provides an numerical solution of Eq |SJ including the 
transient case when the derivative is not equal to zero. 

If the intruder particle moves with a constant velocity the flux is itself (on average) 
constant. A waiting time consistent with EqEJis then generated: 

At = -ln(a)M^), (10) 
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where < £1 < 1 is a uniform random number. 

Given a collision at time t, the probability that the collision occurs on the right is 
given by 

p(+\t) = <f>+(vi)/<t>(vi). (11) 

This is sampled by generating a second uniform random number < £2 < 1- If 
£2 < p{+\t) the collision is on the right hand side: otherwise it takes place on the 
left hand side. 

Having chosen the side, it is then necessary to sample the velocity of the bath 
particle that collides with this side. The probability distribution function of the colliding 
particle's velocity depends on the collision side: 

9+ (v,V 1 )={ fa -«)/(«> if ^l (12) 

otherwise 



g.(v, Vl ) = { (« <0/*-(«i) if«>«i (13) 
otherwise 

The results presented so far apply to any bath velocity distribution - and the 
behavior depends strongly on the exact form 12 . For a Gaussian 



f(v , a) = J — exp(— av 2 ), — 00 < v < 00 (14) 
where a = m/(2kBTB)- The collision fluxes on each side of the tracer particle are 

0±K) = tt(±ui + -i= exp(-ai;j[) +v 1 erf(viy/a)), (15) 

and the total collision rate is 

Mvi) = p(—j= exp(-av 2 ) + vieiiiviyfa)). (16) 
\/7ra 



This function is shown in Figure ^ 

The velocity distribution of the colliding particles, Eq i n the case of a bath 
particle Gaussian velocity distribution is plotted in FigElfor several values of the intruder 
particle velocity. Sampling of this distribution is accomplished with an acceptance- 
rejection method in which the sampling region is adapted to the velocity of the intruder. 
A rare, but possible, case is to select a collision on the right hand side when the surface 
is moving rapidly to the left. The distribution of colliding particles is sharply peaked at 
V\ and it is necessary to take the range of v from v\ — 2 to V\. If the surface is moving 
to the right, the distribution of colliding particle's velocity is more symmetric and one 
can sample v in the range —3 < v < 3. 

Finally, when the intruder particle moving with a velocity v\ collides with a bath 
particle of velocity v the velocity of the former changes instantaneously to 

v 'i = v ^YT^ {v ~ Vxl (17) 
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Figure 1. Flux, or collision rate, on both sides of a surface moving in a bath of particles 
with a Gauss' ' ...... -bottom. 
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Figure 2. Probability distribution of particles that collide with the right hand side of 
a surface moving with a velocity v\ — — 2,-1, 0, 1,2, left to right. The bath particles 
have a Gaussian velocity distribution 

where < a < 1 is the coefficient of restitution and we have taken m = 1 for 
convenience. 

The complete algorithm describing one iteration can now be summarized using 
pseudocode: 

1 Generate a waiting time using EqlTUl 

2 t -> t + At 
while (t out < t) 

{tout tout + 8t 

accumulate averages} 

3 Choose the collision side using Eq ITT1 

4 The velocity of the colliding bath particle v is sampled from the distribution given 
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Figure 3. Ratio of the intruder granular temperature to the bath temperature. 
Symbols show the simulation results and the solid curve shows the theoretical 
prediction of Martin and Piasecki, Eq^] Results are for M — m and n co i = 500000 



by Eq[T21or EqUJ 

5 The post-collisional velocity of the intruder particle is determined from Eq 1171 

Although the time increment between events is variable, averages (mean square 
velocities and velocity distributions) must be computed at equal time intervals. In step 
2, t denotes the total elapsed time and 5t is the constant time interval between the 
accumulation of quantities to be averaged. A convenient choice for St is the average 
collision time. 

Martin and Piasecki obtained an analytic solution of the homogeneous stationary 
Boltzmann equation and showed that the velocity distribution of the intruder particle 
in the steady state is Gaussian and characterized by a temperature, T, that is different 
from the bath temperature Tb- Specifically, the two are related by: 

T _ 1 + a n8 N 
Tb 2 + (l-a)/M' 1 ' 

so that for a < 1, T < Tb- Figure El shows an excellent agreement of the simulation 

with this exact result. 



Since the existence o: 
has been shown recently 



solutions of the Boltzmann equation with power-law tails 



15 



16 



it is interesting to investigate this phenomenon in 
the intruder particle system using the Gillespie algorithm. Therefore, we consider the 
case where the bath particle velocity distribution function takes the following power-law 
form: 



f(a,v) 



2a 



(19) 



7T 1 + a?v A 

The granular temperature of the bath is well defined since the average of the square 
velocity is finite, < v 2 >= l/a. 

Although the Boltzmann equation can no longer be solved in general for an arbitrary 
bath particle velocity distribution (unlike the just-discussed Gaussian distribution) an 
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Figure 4. Granular temperatures of the intruder particle normalized by the 
bath temperature, T/Tb, versus the coefficient of restitution a when the velocity 
distribution of the bath is given by Eq. Ijl9(l for M = 1 and M = 0.5, from top 
to bottom. Dashed curves correspond to the analytical result |2( when the velocity 
distribution of the bath is Gaussian. 



exact solution is possible when the mass ratio is equal to the coefficient of restitution, 
M/m = a. For this specific case one can show that the stationary solution of the 
intruder particle is exactly given by Eq. (JT9j) 3] with a granular temperature equal to the 
bath temperature multiplied by the coefficient of restitution. 

Figure 0] shows the variation of the granular temperature of the intruder particle 
with a for M/m = 1 and M/m = 0.5. When the intruder is light, M < m, the granular 
temperature of the intruder particle is close to the result obtained with the Gaussian 
bath when a > M/m (see Fig. HJ). 

Figure El displays the intruder particle velocity distribution function for different 
values of the coefficient of restitution 0.0 < a < 1.0 when M — 1/2. The exact solution is 
known for a = 0.5, i.e. Eq. (jl9|) with a granular temperature equal to a. The simulation 
results show that in all cases the velocity distribution functions exhibit a power-law tail 
(See inset of FigEJ) with an exponent independent of a and equal to —4. 

3.2. Needle 

In this application a needle intruder, confined to a two-dimensional plane, is immersed in 
a fluid of point particles, each of mass m, at a density p. The needle is characterized by 
its mass M, length L and moment of inertia I and its state is specified by its angular and 
center of mass velocities, ui and Vi, respectively (see FigEJ). The velocity distribution 
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Figure 5. Velocity distribution function of the intruder particle for different values of 
the coefficient of restitution a — 0.0,0.1,0.2, ...0.8,0.9, 1.0, from top to bottom, when 
M = 1/2 and when the bath distribution is given by Ea. H19|) . The insert is a log- log 
plot showing the algebraic decay of the velocity distribution 




Figure 6. Geometry of the needle-point system 



of the point particles is again given by Eqll4l 

Two main modifications of the Gillespie algorithm are required in order to simulate 
this system. First, if v\ ^ the flux is not (on average) constant between collisions. 
Second, it is necessary to select the point of impact on the needle. 

The collision flux on both sides of the needle at a point —L/2 < x < L/2 is 
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0(vi.n + uox). Unlike the case considered above, this flux is time-dependent if u 7^ 
since the normal vector n rotates. Specifically 

v n = Vi.n(t) = vi y cos{ut + 9o) — v\ x sin(ut + 9 ), (20) 

where 8q is the orientation of the needle at t = 0. The total flux over the entire length 
of the needle is 

rL/2 

$(£) = / <f)(w 1 .n(t) +ux)dx. (21) 

J-L/2 

For the sake of simplicity we take a = 1 and L = 1 and we obtain 



f vl u v n 1 \ „ / u\ 



e n i 
+ P- 



/ i _ 2t>n ) _|_ e -^ w ( 1 + — 



A collision time is generated by solving the equation 



/ df$(f)=-]n(£3 
Jo 



(22) 



(23) 



Since the integral cannot be performed analytically, we use the Newton-Raphson 
method: 

J^(t')dt' + H^) 

ti ~ t W) — ( } 



with an initial guess t = — ^ 3 ) • The procedure is iterative, i.e. t is substituted by t\ in 

Eq.(j24j), etc. until | /J 1 (^(t')dt' +\n(^)\ is smaller than the required precision. In general 
the convergence is fast and only a few iterations are required. Occasionally, when $ is 
small, the Newton-Raphson method oscillates between two "stable" positions and there 
is no convergence. When this situation arises, we switch to a bisection procedure. With 
this modification, the method seems to be robust. 

Once the time to collision has been selected, one updates the normal velocity of the 
center-of-mass of the needle using Eq EH 

In order to choose the position of the impact, one has to calculate the probability 
that the collision occurs at a distance x from the center of mass of the needle, whatever 
the velocities of the bath particles, at a given time t 

P(x\t)= L" -1/2 < x < 1/2 (25) 

This probability is clearly not uniform over the length of the needle. But since it is 
concave, the maximum is obtained for x = 1/2 or x = —1/2 (if v n = the probability is 
a maximum at both ends). We select x using a standard acceptance- rejection method. 
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Once the position of the impact has been decided, one has to choose if the bath 
particle collides on the right or left hand side of the needle. The probability of the 
former event is 

4w + «*) (26) 

One chooses a random number between and 1. If £ < p(+\x,t), the particle collides 
with a particle from the right, otherwise the collision is on the left. 

Next the velocity of the colliding bath particle must be selected. The probability 
that the colliding particle has a velocity between v and v + dv is g±(v, v n + xuS)dv. It is 
even more important to sample this distribution carefully than in the one-dimensional 
case as more extreme velocities are encountered in this system. 

Finally, the needle velocity and angular velocity are updated using 

vi = vi + n— (27) 

Iuj' = Iu + xAp (28) 

where 

(1 + a)g.n . , 

m ^ M ' I 

n is a unit vector normal to the length of the needle and g = vi — v + cuxn is the relative 
velocity at the point of impact (the velocity of the colliding bath particle also changes, 
but we do not need to know the new value). 

It is necessary to simulate a few hundred thousand to several million collisions in 
order to obtain good estimates for the properties of interest. The convergence of the 
simulation is slower when the mass of the needle is much bigger than the bath particle 
mass. 

We have previously used this algorithm to confirm a kinetic theory prediction that, 
when the coefficient of restitution is smaller than unity, the temperature of the bath 
is larger than the translational gran ular temperature which is in turn larger than the 
rotational granular temperature [lQj . 

Here we present new results for the cross-correlation between the two degrees of 
freedom of the needle as a function of the coefficient of restitution, a, and for different 
values of the mass ratio Mj m. Specifically, we have calculated 

< v 2 uj 2 > , , 

R = 5 5—, 30 

< V 2 >< UJ 2 > 

which is equal to one in equilibrium systems and obviously independent of the mass 
ratio. Conversely, one observes that for small mass ratios in an inelastic system, there 
is a positive correlation that increases with decreasing a: See Figure We expect this 
to be a generic feature of anisotropic granular particles in any dimension, regardless of 
the bath particle velocity distribution. 

Unlike the one- dimensional model, where the velocity distribution of the intruder 
particle is Gaussian for all values of the coefficient of restitution, the translational and 
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angular velocity distributions of the granular needle are never strictly Gaussian (except 
for elastic collisions). For a large range of values of M/m and a, however, the Gaussian is 
a very accurate approximation. It is only for a light needle and highly inelastic collisions 
that deviations start to become apparent: See Figure |H1 
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4. DMSC versus Gillepsie algorithm 

The two simulation methods provide an exact numerical solution o 



a Boltzmann 

equation (they can also be used for inelastic Maxwell models 15; 16; El|)- The DSMC 
algorithm has been described in several articles 4; 18], For convenience we describe here 
the version that would be applied to the one- dimensional system discussed in section 
3.1. 

Velocities sampled from a Gaussian distribution, EqH3J are assigned to nbath bath 
particles. A side of the intruder particle is then selected at random: an = +1 for a 
collision on the left, an = —1 for a collision on the right. A bath particle, i, is then 
selected randomly. The collision is accepted with probability Q{vn(Jn)unl ^max where 
6(x) is the Heaviside function, ujh = 2p\vn\ and uj max is an upper bound estimate 
of the probability that a particle collides per unit time. If the collision is accepted a 
post-collisional velocity, computed from Eq[T7l is assigned to the intruder particle. If 
ujh > oj max the estimate of the latter is updated: u max <— uJn- 

We have implemented this algorithm for the one-dimensional intruder described in 
section 3.1. We obtain the same results with comparable computational effort. 

5. Conclusion 

We have shown how the Gillespie algorithm can be used to obtain an exact numerical 
solution of the Boltzmann equation of an intruder particle in a bath of particles with 
an arbitrary velocity distribution. We used the method to obtain new results for 
a one-dimensional system consisting of an intruder particle in a bath with a power 
law distribution. We also used it to demonstrate, for the first time, the presence of 
correlations between the translational and rotational momenta of an anisotropic particle. 

Although the results presented here apply to the steady state, the method is equally 
valid for the transient case. It is clear that the Gillespie algorithm offers no significant 
computational advantage over the DSMC method for these intruder particle systems. 
It is of interest, however, that two apparently dissimilar methods can be applied to the 
same physical system. Finally, we note that, as with DSMC, the Gillespie method can 
be easily generalized to three-dimensional systems. 
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